C***************************************************************
C             Subroutine LSEMA by Stephen Kirkup                     
C***************************************************************
C 
C  Copyright 2001- Stephen Kirkup
C  School of Computing Engineering and Physical Sciences
C  University of Central Lancashire - www.uclan.ac.uk 
C  smkirkup@uclan.ac.uk
C  http://www.researchgate.net/profile/Stephen_Kirkup
C
C  This open source code can be found at
C   www.boundary-element-method.com/fortran/LSEMA.FOR 
C
C  Issued under the GNU General Public License 2007, see gpl.txt

C
C  Part of the the author's open source BEM packages. 
C  All codes and manuals can be downloaded from 
C  www.boundary-element-method.com
C
C***************************************************************
C This subroutine computes the solution of Laplace's equation 
C                  __ 2            
C                  \/    {\phi}     =  0   
C
C  in the region exterior to a set of thin surfaces.
C
C
C How to use the subroutine
C -------------------------
C
C The following diagram shows how the subroutine is to be used. 
C
C                                   ....................................
C                                   :                                  :
C                                   :                                  :
C      ----------------------       :     --------------------------   :
C      |                    |       :     |                        |   :
C      |   MAIN PROGRAM     |------->-----|         LSEMA          |   :
C      | (e.g. LSEMA_t.for) |       :     |                        |   :
C      |                    |       :     --------------------------   :
C      ----------------------       :                 |                :
C                                   :                 >                :
C                                   :                 |                :
C                                   :      ------------------------    :
C          Package ---------------->:      | subordinate routines |    :
C                                   :      ------------------------    :
C                                   :                                  :
C                                   :      (this file)                 :  
C                                   :..................................:
C                                  /         |                 |
C                               |_           >                 >
C                              /             |                 |
C             ................       ................   ................  
C             :              :       :   --------   :   :  ---------   : 
C             : (geom3d.for) :---<---:   |L3ALC |   :   :  | GLS   |  : 
C             : (geom2d.for) :       :   --------   :   :  ---------   :  
C             :..............:       : -------------:   : -------------:  
C                                    : |subordinate|:   : |subordinate|: 
C                                    : | routines  |:   : | routines  |:  
C                                    : -------------:   : -------------: 
C                                    :              :   :              : 
C                                    : (l3alc.for)  :   :  (gls.for)   :
C                                    :..............:   :..............:
C                                    
C
C The contents of the main program must be linked to LSEMA.FOR, L3ALC.FOR
C GLS2.FOR, GL8.FOR, GEOM3D.FOR and GEOM2D.FOR.
C
C Method of solution
C ------------------
C 
C In the main program, the shell(s) must be described as a set of
C  panels. The panels are defined by two indices (integers) which
C  label a node or vertex on the shell generator. The data structure
C  VERTEX lists and enumerates the coordinates of the vertices, the data
C  structure HELV defines each panel by indicating the labels for
C  the three nodes that are its vertices and hence enumerates the
C  panels.
C                                 
C Format and parameters of the subroutine
C ---------------------------------------
C
C The subroutine has the form
C
C      SUBROUTINE LSEMA(MAXNV,NV,VERTEX,MAXNH,NH,HELV,
C     *                 MAXNPE,NPE,PEXT,
C     *                 HA,HB,HF,HAA,HBB,HFF,
C     *                 HIPHI,HIVEL,PEIPHI,
C     *                 LSOL,LVALID,EGEOM,
C     *                 PHIDIF,PHIAV,VELDIF,VELAV,PEPHI,
C     *                 AMAT,BMAT,L_EH, M_EH,
C     *                 PERM, XORY,C,WKSPC1,WKSPC2,WKSPC3)

C The parameters to the subroutine
C ================================

C Geometry of the shell {\Pi} (input)
C integer MAXNV: The limit on the number of vertices of the triangles
C  that defines (approximates) {\Pi}. MAXNV>=2.
C integer NV: The number of vertices on {\Pi}. 2<=NV<=MAXNV.
C real VERTEX(MAXNV,2): The coordinates of the vertices. VERTEX(i,1),
C  VERTEX(i,2) are the x,y,z coordinates of the i-th vertex. 
C integer MAXNH: The limit on the number of panels describing {\Pi}.
C  MAXNH>=1.
C integer NH: The number of panels describing {\Pi}. 1<=NH<=MAXNH.
C integer HELV(MAXNH,2): The indices of the two vertices on the
C  generator defining each panel. The i-th panel has vertices 
C  (VERTEX(HELV(i,1),1),VERTEX(HELV(i,1),2))) and
C  (VERTEX(HELV(i,2),1),VERTEX(HELV(i,2),2))).

C Exterior points at which the solution is to be observed (input)
C integer MAXNPE: Limit on the number of points exterior to the
C  shell. MAXNPE>=1.
C integer NPE: The number of exterior points. 0<=NPE<=MAXNPE.
C real PEXT(MAXNPE,2). The coordinates of the exterior point.
C  PEXT(i,1),PEXT(i,2) are the x,y,z coordinates of the i-th point. 

C Shell Conditions
C  a(p){\delta}(p)+b(p){\nu}(p)=f(p)
C  real HA(MAXNH).  Function a at the central points
C  real HB(MAXNH).  Function b at the central points
C  real HF(MAXNH).  Function f at the central points
C  A(p){\Phi}(p)+B(p)V(p)=F(p)
C  real HAA(MAXNH).  Function A at the central points
C  real HBB(MAXNH).  Function B at the central points
C  real HFF(MAXNH).  Function F at the central points

C Incident field
C real HIPHI(MAXNH). The incident potential {\phi} at the central
C  points
C real HIVEL(MAXNH). The derivative of the incident potential {\phi} 
C  at the  central points
C real PEIPHI(MAXNPE). The incident potential {\phi} at the exterior 
C  points

                  
C Validation and control parameters (input) 
C logical LSOL: A switch to control whether the particular solution is
C  required.
C logical LVALID: A switch to enable the choice of checking of 
C  subroutine parameters.
C real EGEOM: The maximum absolute error in the parameters that
C  describe the geometry.

C Solution (output)
C real SPHI(MAXNH): The potential ({\phi}) at the centres of the 
C  boundary panels.
C real SVEL(MAXNH): The (v or d{\phi}/dn where n is the outward 
C  normal to the boundary) at the centres of the boundary 
C  panels.
C real PEPHI(MAXNPE): The potential ({\phi}) at the exterior points.

C  Working space 
C      REAL*8     AMAT(2*MAXNH,2*MAXNH)
C      REAL*8     BMAT(2*MAXNH,2*MAXNH)
C      REAL*8     L_EH(MAXNPE,MAXNH)
C      REAL*8     M_EH(MAXNPE,MAXNH)
C      REAL*8     C(2*MAXNH)
C      REAL*8     HALPHA(2*MAXNH)
C      REAL*8     HBETA(2*MAXNH)
C      REAL*8     HFFF(2*MAXNH)
C      REAL*8     PHIINF(2*MAXNH)
C      REAL*8     VELINF(2*MAXNH)
C      REAL*8     WKSPC1(2*MAXNH)
C      REAL*8     WKSPC2(2*MAXNH)
C      REAL*8     WKSPC3(2*MAXNH)

C Notes on the geometric parameters
C ---------------------------------
C (1) Each of the vertices listed in VERTEX must be distinct points
C  with respect to EGEOM.
C (2) The boundary must be complete and closed. Thus 
C  HELV(i,2)=HELV(i+1,1) for i=1..NH-1 and VERTEX(HELV(1,1),1)=0
C  and VERTEX(HELV(NH,2),1)=0.
C (3) The indices of the nodes listed in HELV must be such that they
C  are ordered counter-clockwise around the generator of the boundary.
C (4) The generator of the largest panel must be no more than 10x 
C  the length of the generator of the smallest panel.

C Notes on the exterior points 
C ----------------------------
C (1) The points in PEXT should lie outside the shell, as defined
C  by the parameters VERTEX and HELV. Any point lying outside the 
C  shell.


C The subroutine

      SUBROUTINE LSEMA(MAXNV,NV,VERTEX,MAXNH,NH,HELV,
     *                 MAXNPE,NPE,PEXT,
     *                 HA,HB,HF,HAA,HBB,HFF,
     *                 HIPHI,HIVEL,PEIPHI,
     *                 LSOL,LVALID,EGEOM,
     *                 PHIDIF,PHIAV,VELDIF,VELAV,PEPHI,
     *                 AMAT,BMAT,L_EH, M_EH,
     *                 PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3)
      PARAMETER (MAXNGQ=32)
      PARAMETER (MAXNTQ=10000)

C  Boundary geometry
C   Limit on the number of vertices on {\Pi}
      INTEGER    MAXNV
C   The number of vertices on {\Pi}
      INTEGER    NV
C   The coordinates of the vertices on {\Pi}
      REAL*8     VERTEX(MAXNV,2)
C   Limit on the number of panels describing {\Pi}
      INTEGER    MAXNH
C   The number of panels describing {\Pi}
      INTEGER    NH
C   The indices of the vertices describing each panel
      INTEGER    HELV(MAXNH,2)
      
C  Exterior points at which the solution is to be observed
C   Limit on the number of points exterior to the shell where 
C    solution is sought
      INTEGER    MAXNPE
C   The number of exterior points
      INTEGER    NPE
C   Coordinates of the exterior points
      REAL*8     PEXT(MAXNPE,2)


C  Shell Conditions
C   a(p){\delta}(p)+b(p){\nu}(p)=f(p)
C    function a at the central points
      REAL*8     HA(MAXNH)
C    function b at the central points
      REAL*8     HB(MAXNH)
C    function f at the central points
      REAL*8     HF(MAXNH)
C   A(p){\Phi}(p)+B(p)V(p)=F(p)
C    function A at the central points
      REAL*8     HAA(MAXNH)
C    function B at the central points
      REAL*8     HBB(MAXNH)
C    function F at the central points
      REAL*8     HFF(MAXNH)

C  Incident field
C    The incident potential {\phi} at the central points
      REAL*8    HIPHI(MAXNH)
C    The derivative of the incident potential {\phi} at the 
C     central points
      REAL*8    HIVEL(MAXNH)
C    The incident potential {\phi} at the exterior points
      REAL*8    PEIPHI(MAXNPE)

C Validation and control parameters
      LOGICAL   LSOL
      LOGICAL   LVALID
      REAL*8    EGEOM

C Solution
      REAL*8    PHIDIF(MAXNH)
      REAL*8    PHIAV(MAXNH)
      REAL*8    VELDIF(MAXNH)
      REAL*8    VELAV(MAXNH)
      REAL*8    PEPHI(MAXNPE)

C  Working space 
      REAL*8    AMAT(2*MAXNH,2*MAXNH)
      REAL*8    BMAT(2*MAXNH,2*MAXNH)
      REAL*8    L_EH(MAXNPE,MAXNH)
      REAL*8    M_EH(MAXNPE,MAXNH)
      REAL*8    C(2*MAXNH)
      REAL*8    HALPHA(2*MAXNH)
      REAL*8    HBETA(2*MAXNH)
      REAL*8    HFFF(2*MAXNH)
      REAL*8    PHIINF(2*MAXNH)
      REAL*8    VELINF(2*MAXNH)
      REAL*8    WKSPC1(2*MAXNH)
      REAL*8    WKSPC2(2*MAXNH)
      REAL*8    WKSPC3(2*MAXNH)

C Further information from system solution GLS
      INTEGER*4 PERM(2*MAXNH)
      LOGICAL   XORY(2*MAXNH)

      REAL*8    DIST2
      REAL*8    DIAM

C  Constants
C   Real scalars: 0, 1, 2, half, pi
      REAL*8    ZERO,ONE,TWO,PI


C  Geometrical description of the shell
C   panels counter
      INTEGER   IPI,JPI
C   The points exterior to the shell where the solution is sought 
      INTEGER   IPE
C   Parameters for L3ALC
      REAL*8    P(2),PA(2),PB(2),QA(2),QB(2),QC(2),VECP(2)
      LOGICAL   LPONEL

C  Quadrature rule information
C   [Note that in this program two quadrature rules are used: one for
C    the case when the point P lies on the panel (LPONEL=.TRUE.) and
C    one for the case when P does not lie on the panel.]
C   Quadrature rule used when LPONEL=.TRUE.
C    Number of quadrature points
      INTEGER   NGQON
C    Abscissae of the actual quadrature rule
      REAL*8    AGQON(MAXNGQ)
C    Weights of the actual quadrature rule
      REAL*8    WGQON(MAXNGQ)
C   Quadrature rule used when LPONEL=.FALSE.
C    Number of quadrature points
      INTEGER   NGQOFF
C    Abscissae of the actual quadrature rule
      REAL*8    AGQOFF(MAXNGQ)
C    Weights of the actual quadrature rule
      REAL*8    WGQOFF(MAXNGQ)
C   Quadrature rule parameters for L3ALC
C    Actual number of quadrature points
      INTEGER   NGQ
C    Abscissae of the actual quadrature rule in the generator direction
      REAL*8    AGQ(MAXNGQ)
C    Weights of the actual quadrature rule
      REAL*8    WGQ(MAXNGQ)
C   Counter through the quadrature points
      INTEGER   IGQ
C    Abscissae of the actual quadrature rule in the theta direction
      REAL*8    ATQ(MAXNTQ)
C    Weights of the actual quadrature rule
      REAL*8    WTQ(MAXNTQ)


C  Validation and control parameters for subroutine L3ALC
      LOGICAL   LVAL
      REAL*8    EQRULE
      LOGICAL   LL
      LOGICAL   LM
      LOGICAL   LMT
      LOGICAL   LN

C  Parameters for subroutine L3ALC. 
      REAL*8    DISL
      REAL*8    DISM
      REAL*8    DISMT
      REAL*8    DISN

C  Other variables
C   Failure flag
      LOGICAL   LFAIL
C   Accumulation of solution {\phi}
      REAL*8    SUMPHI
C   Error flag
      LOGICAL   LERROR     

      REAL*8    WKSPCE(2*MAXNTQ+MAXNGQ)
      REAL*8    SIZMAX,SIZMIN

C INITIALISATION
C --------------

C Set constants
      ZERO=0.0D0
      ONE=1.0D0
      TWO=2.0D0
      PI=4.0D0*ATAN(ONE)


C Set up validation and control parameters
C  Switch off the validation of L3ALC
      LVAL=.FALSE.
C  Set EQRULE
      EQRULE=1.0D-6
C  Set EGEOM
      EGEOM=1.0D-6

C Set up the quadrature rule(s).
C  Set up quadrature rule for the case when P is not on the panel.
C   Set up 8 point Gauss-Legendre rules
      CALL GL8(MAXNGQ,NGQOFF,WGQOFF,AGQOFF)
C  Set up quadrature rule for the case when P is on the panel.
C   This is done by splitting the quadrature rule at the centre.
      NGQON=2*NGQOFF
      DO 330 IGQ=1,NGQOFF
        AGQON(IGQ)=AGQOFF(IGQ)/TWO
        AGQON(NGQOFF+IGQ)=0.5D0+AGQOFF(IGQ)/TWO
        WGQON(IGQ)=WGQOFF(IGQ)/TWO
        WGQON(NGQOFF+IGQ)=WGQOFF(IGQ)/TWO
330   CONTINUE



C Validation
C ==========

C Validation of parameters of LEBEM3
C ---------------------------------

      IF (LVALID) THEN

C Validation of main paramters
        LERROR=.FALSE.
        IF (MAXNV.LT.2) THEN
          WRITE(*,*) 'MAXNV = ',MAXNV
          WRITE(*,*) 'ERROR(LSEMA) - must have MAXNV>=4'
          LERROR=.TRUE.
        END IF
        IF (NV.LT.2.OR.NV.GT.MAXNV) THEN
          WRITE(*,*) 'NV = ',NV
          WRITE(*,*) 'ERROR(LSEMA) - must have 4<=NV<=MAXNV'
          LERROR=.TRUE.
        END IF
        IF (MAXNH.LT.2) THEN
          WRITE(*,*) 'MAXNH = ',MAXNH
          WRITE(*,*) 'ERROR(LSEMA) - must have MAXNH>=2'
          LERROR=.TRUE.
        END IF
        IF (NH.LT.3.OR.NH.GT.MAXNH) THEN
          WRITE(*,*) 'NH = ',NH
          WRITE(*,*) 'ERROR(LSEMA) - must have 4<=NH<=MAXNH'
          LERROR=.TRUE.
        END IF
        IF (MAXNPE.LT.1) THEN
          WRITE(*,*) 'MAXNPE = ',MAXNPE
          WRITE(*,*) 'ERROR(LSEMA) - must have MAXNPE>=1'
          LERROR=.TRUE.
        END IF
        IF (NPE.LT.0.OR.NPE.GT.MAXNPE) THEN
          WRITE(*,*) 'NPE = ',NPE
          WRITE(*,*) 'ERROR(LSEMA) - must have 3<=NPE<=MAXNPE'
          LERROR=.TRUE.
        END IF
        IF (EGEOM.LE.ZERO) THEN
          WRITE(*,*) 'NPE = ',NPE
          WRITE(*,*) 'ERROR(LSEMA) - EGEOM must be positive'
          LERROR=.TRUE.
        END IF
        IF (LERROR) THEN
          LFAIL=.TRUE.
          WRITE(*,*)
          WRITE(*,*) 'Error(s) found in the main parameters of LSEMA'
          WRITE(*,*) 'Execution terminated'
          STOP
        END IF
      END IF

C Find the diameter DIAM of the boundary
      DIAM=0.0
      DO 100 IV=1,NV-1
        PA(1)=VERTEX(IV,1)
        PA(2)=VERTEX(IV,2)
        DO 110 JV=IV+1,NV
          PB(1)=VERTEX(JV,1)
          PB(2)=VERTEX(JV,2)
          DIAM=MAX(DIAM,DIST2(PA,PB))
110     CONTINUE
100   CONTINUE

      IF (LVALID) THEN
        LERROR=.FALSE.
C Check that the boundary defined by HELV is complete and closed
        DO 120 IPI=1,NH
120     CONTINUE


C Check that EGEOM is not too large
        IF (EGEOM.GT.DIAM/100.0D0) THEN
          WRITE(*,*) 'EGEOM = ',EGEOM
          WRITE(*,*) 'ERROR(LSEMA) - EGEOM is set too large'
          LERROR=.TRUE.
        END IF
        IF (LERROR) THEN
          LFAIL=.TRUE.
          WRITE(*,*)
          WRITE(*,*) 'Error in boundary geometry or EGEOM'
          WRITE(*,*) 'Execution terminated'
        END IF
      END IF                  

      IF (LVALID) THEN
C Check that the vertices are distinct with respect to EGEOM
        LERROR=.FALSE.
        DO 130 IV=1,NV-1
          PA(1)=VERTEX(IV,1)
          PA(2)=VERTEX(IV,2)
          DO 140 JV=IV+1,NV
            PB(1)=VERTEX(JV,1)
            PB(2)=VERTEX(JV,2)
            IF (ABS(PA(1)-PB(1)).LT.EGEOM) THEN
              IF (ABS(PA(2)-PB(2)).LT.EGEOM) THEN
                WRITE(*,*) 'Vertices ',IV,JV,' are not distinct'
                LERROR=.TRUE.
              END IF
            END IF
140       CONTINUE
130     CONTINUE
        IF (LERROR) THEN
          WRITE(*,*) 
          WRITE(*,*) 'ERROR(LSEMA) - Vertices (see above) coincide'
          WRITE(*,*) 'Execution terminated'
          STOP
        END IF
      END IF          


C Check that the panels are not of disproportionate sizes
      IF (LVALID) THEN
        SIZMAX=ZERO
        SIZMIN=DIAM
        DO 150 IPI=1,NH
          QA(1)=VERTEX(HELV(IPI,1),1)
          QA(2)=VERTEX(HELV(IPI,1),2)
          QB(1)=VERTEX(HELV(IPI,2),1)
          QB(2)=VERTEX(HELV(IPI,2),2)
          SIZMAX=MAX(SIZMAX,DIST2(QA,QB))
          SIZMIN=MIN(SIZMIN,DIST2(QA,QB))
150     CONTINUE
        IF (SIZMAX.GT.10.0D0*SIZMIN) THEN
          WRITE(*,*) 'WARNING(LSEMA) - panels of disproportionate'
          WRITE(*,*) ' sizes'
        END IF
      END IF          
    
          
C Check that the boundary does not contain sharp angles
      IF (LVALID) THEN
        LERROR=.FALSE.
        DO 160 IPI=2,NH
          IF (HELV(IPI-1,2).EQ.HELV(IPI,1)) THEN
            QA(1)=VERTEX(HELV(IPI,1),1)
            QA(2)=VERTEX(HELV(IPI,1),2)
            QB(1)=VERTEX(HELV(IPI,2),1)
            QB(2)=VERTEX(HELV(IPI,2),2)
            QC(1)=VERTEX(HELV(IPI-1,1),1)
            QC(2)=VERTEX(HELV(IPI-1,1),2)
            IF (DIST2(QC,QB).LT.MAX(DIST2(QA,QB),DIST2(QA,QC))) THEN
              WRITE(*,*) 'Sharp angle at node ',HELV(IPI,1)
              LERROR=.TRUE.
            END IF
          END IF
160     CONTINUE
        IF (LERROR) THEN
          WRITE(*,*)
          WRITE(*,*) 'WARNING(LSEMA) - Boundary has sharp angles'
        END IF
      END IF          
     

C Validation of the surface functions
      IF (LVALID.AND.LSOL) THEN
        LERROR=.FALSE.
        DO 170 IPI=1,NH
          IF (MAX(ABS(HA(IPI)),ABS(HB(IPI))).LT.1.0D-6) 
     *     LERROR=.TRUE.
170     CONTINUE
        IF (LERROR) THEN
          WRITE(*,*) 
          WRITE(*,*) 'ERROR(LSEMA) - at most one of HA(i),HB(i)'
          WRITE(*,*) ' may be zero for all i'
          WRITE(*,*) 'Execution terminated'
          STOP
        END IF
      END IF
        
C  Compute the discrete L, M, Mt and N matrices
C   Loop(IPI) through the points on the boundary
      DO 510 IPI=1,NH
C    Set P
        PA(1)=VERTEX(HELV(IPI,1),1)
        PA(2)=VERTEX(HELV(IPI,1),2)
        PB(1)=VERTEX(HELV(IPI,2),1)
        PB(2)=VERTEX(HELV(IPI,2),2)
        P(1)=(PA(1)+PB(1))/TWO
        P(2)=(PA(2)+PB(2))/TWO
C    Set VECP to the normal on the boundary of the panel at P
        CALL NORM2(PA,PB,VECP)
C    Loop(IPI) through the panels
        DO 520 JPI=1,NH
C     Set QA and QB, the coordinates of the edges of the JSEth panel
          QA(1)=VERTEX(HELV(JPI,1),1)
          QA(2)=VERTEX(HELV(JPI,1),2)
          QB(1)=VERTEX(HELV(JPI,2),1)
          QB(2)=VERTEX(HELV(JPI,2),2)

C     Set LPONEL
          IF (IPI.EQ.JPI) THEN
            LPONEL=.TRUE.
          ELSE
            LPONEL=.FALSE.
          END IF

C     Select quadrature rule for L3ALC
C   :  Select the quadrature rule AGQON-WGQON in the case when the
C   :   point p lies on the panel, otherwise select AGQOFF-WGQOFF
C      [Note that the overall method would benefit from selecting from
C       a wider set of quadrature rules, and an appropriate method
C       of selection]
          IF (LPONEL) THEN
          NGQ=NGQON
          DO 600 IGQ=1,NGQ
            AGQ(IGQ)=AGQON(IGQ)
            WGQ(IGQ)=WGQON(IGQ)
600       CONTINUE
          ELSE 
          NGQ=NGQOFF
          DO 610 IGQ=1,NGQ
            AGQ(IGQ)=AGQOFF(IGQ)
            WGQ(IGQ)=WGQOFF(IGQ)
610       CONTINUE
          END IF



C Quadrature rule in the theta direction is constructed out of individual
C Gauss rules so that the length of each is approximately equal to the
C length of the panel at the generator.
          RADMID=(QA(1)+QB(1))/TWO
          SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+
     *     (QA(2)-QB(2))*(QA(2)-QB(2))
          GLEN=SQRT(SGLEN)
          CIRMID=PI*RADMID
          NDIV=1+CIRMID/GLEN
          TDIV=ONE/DBLE(NDIV)
          NTQ=NDIV*NGQ
          IF (NTQ.GT.MAXNTQ) THEN
            WRITE(*,*) 'ERROR(LSEMA) - MAXNTQ is set too small'
            STOP
          END IF
          DO 145 IDIV=1,NDIV
            DO 155 IGQ=1,NGQ
              WTQ((IDIV-1)*NGQ+IGQ)=WGQ(IGQ)/DBLE(NDIV)
              ATQ((IDIV-1)*NGQ+IGQ)=AGQ(IGQ)/DBLE(NDIV)+
     *         TDIV*DBLE(IDIV-1)
155         CONTINUE
145       CONTINUE


C    All operators are required
          LL=.TRUE.
          LM=.TRUE.
          LMT=.TRUE.
          LN=.TRUE.


C    Call of L3ALC routine to compute [L], [M], [Mt], [N]
          CALL L3ALC(P,VECP,QA,QB,LPONEL,
     *     MAXNGQ,NGQ,AGQ,WGQ,MAXNTQ,NTQ,ATQ,WTQ,
     *     LVAL,EGEOM,EQRULE,LFAIL,
     *     LL,LM,LMT,LN,DISL,DISM,DISMT,DISN,
     *     WKSPCE)


          AMAT(IPI,JPI)=-DISM
          AMAT(IPI,NH+JPI)=0.0D0
          AMAT(NH+IPI,JPI)=-DISN
          AMAT(NH+IPI,NH+JPI)=0.0D0
          BMAT(IPI,JPI)=-DISL
          BMAT(IPI,NH+JPI)=0.0D0
          BMAT(NH+IPI,JPI)=-DISMT
          BMAT(NH+IPI,NH+JPI)=0.0D0
         
C    Close loop(JPI) 
520     CONTINUE

        AMAT(IPI,NH+IPI)=1.0D0
        BMAT(NH+IPI,NH+IPI)=-1.0D0
            
C   Close loop(IPI) 
510   CONTINUE

C Set incident field to zero
	  DO 550 IPI=1,NH
        HIPHI(IPI)=0.0D0
        HIVEL(IPI)=0.0D0
550   CONTINUE
      DO 560 IPI=1,NPE
        PEIPHI(IPI)=0.0D0
560   CONTINUE


C  Set C,HALPHA,HBETA,HFFF
      DO 700 IPI=1,NH      
C  Set C as the incident potential 
        C(IPI)=HIPHI(IPI)
        C(NH+IPI)=HIVEL(IPI)        
C  Set HALPHA as the coefficient of [{\delta} {\Phi}]
        HALPHA(IPI)=HA(IPI)
        HALPHA(NH+IPI)=HAA(IPI)
C  Set HBETA as the coefficient of [{\nu} V]
        HBETA(IPI)=HB(IPI)
        HBETA(NH+IPI)=HBB(IPI)
C  Set HFFF as [f F]
        HFFF(IPI)=HF(IPI)
        HFFF(NH+IPI)=HFF(IPI)
700   CONTINUE      

C  On solution C becomes {\zeta}/{\epsilon} at the central points
      CALL GLS(2*MAXNH,2*NH,AMAT,BMAT,C,HALPHA,HBETA,
     *          HFFF,PHIINF,VELINF,LFAIL,PERM,XORY,WKSPC1,WKSPC2,WKSPC3)

C  Set C
      DO 710 IPI=1,NH
        PHIDIF(IPI)=PHIINF(IPI)
        PHIAV(IPI)=PHIINF(NH+IPI)
        VELDIF(IPI)=VELINF(IPI)
        VELAV(IPI)=VELINF(NH+IPI)
710   CONTINUE
        
C   Compute sound pressures at the selected exterior points.
C    Loop through the the points in the exterior region
      DO 800 IPE=1,NPE
C    Set P
        P(1)=PEXT(IPE,1)
        P(2)=PEXT(IPE,2)
C    Set VECP, this is arbitrary as the derivative of the potential at P
C     is not sought.
        VECP(1)=ONE
        VECP(2)=ZERO

C    Initialise SUMPHI to the incident potential
        SUMPHI=PEIPHI(IPE)

C    Loop(JPI) through the panels
        DO 850 JPI=1,NH
C     Compute the discrete L and M integral operators. 
            
C     Set QA and QB, the coordinates of the edges of the JSEth panel
          QA(1)=VERTEX(HELV(JPI,1),1)
          QA(2)=VERTEX(HELV(JPI,1),2)
          QB(1)=VERTEX(HELV(JPI,2),1)
          QB(2)=VERTEX(HELV(JPI,2),2)

C Quadrature rule in the generator direction
          NGQ=NGQOFF
          DO 1610 IGQ=1,NGQ
            AGQ(IGQ)=AGQOFF(IGQ)
            WGQ(IGQ)=WGQOFF(IGQ)
1610      CONTINUE

C Quadrature rule in the theta direction is constructed out of individual
C Gauss rules so that the length of each is approximately equal to the
C length of the panel at the generator.
          RADMID=(QA(1)+QB(1))/TWO
          SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+
     *     (QA(2)-QB(2))*(QA(2)-QB(2))
          GLEN=SQRT(SGLEN)
          CIRMID=PI*RADMID
          NDIV=1+CIRMID/GLEN
          TDIV=ONE/DBLE(NDIV)
          NTQ=NDIV*NGQ
          IF (NTQ.GT.MAXNTQ) THEN
            WRITE(*,*) 'ERROR(LSEMA) - MAXNTQ is set too small'
            STOP
          END IF
          DO 1145 IDIV=1,NDIV
            DO 1155 IGQ=1,NGQ
              WTQ((IDIV-1)*NGQ+IGQ)=WGQ(IGQ)/DBLE(NDIV)
              ATQ((IDIV-1)*NGQ+IGQ)=AGQ(IGQ)/DBLE(NDIV)+
     *         TDIV*DBLE(IDIV-1)
1155        CONTINUE
1145      CONTINUE

C     All the points do not lie on the boundary hence LPONEL=.FALSE.
          LPONEL=.FALSE.              

C     Only L, M operators are required. Set LL,LM true, 
C      LMT,LN false. 
          LL=.TRUE.
          LM=.TRUE.
          LMT=.FALSE.
          LN=.FALSE.
                
C     Call L3ALC.
             CALL L3ALC(P,VECP,QA,QB,LPONEL,
     *        MAXNGQ,NGQ,AGQ,WGQ,MAXNTQ,NTQ,ATQ,WTQ,
     *        LVAL,EGEOM,EQRULE,LFAIL,
     *        LL,LM,LMT,LN,DISL,DISM,DISMT,DISN,
     *        WKSPCE)


C     Accumulate phi 
          SUMPHI=SUMPHI+DISM*PHIDIF(JPI)-DISL*VELDIF(JPI)
          L_EH(IPE,JPI)=DISL
          M_EH(IPE,JPI)=DISM

C      Close loop (JPI) through the panels
850     CONTINUE

        PEPHI(IPE)=SUMPHI

C     Close loop(IPE) through the exterior points
800   CONTINUE

      END

C ----------------------------------------------------------------------

C Subordinate routines for LEBEMA
C ===============================

C Subroutines required for L3ALC (not in file L3ALC.FOR) 
C  Subroutine for returning the square root.
       REAL*8 FUNCTION FNSQRT(X)
       REAL*8 X
       FNSQRT=SQRT(X)
       END